library(tidyverse)
Registered S3 methods overwritten by 'dbplyr':
method from
print.tbl_lazy
print.tbl_sql
── Attaching packages ────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.3 ✓ purrr 0.3.4
✓ tibble 3.1.2 ✓ dplyr 1.0.6
✓ tidyr 1.1.3 ✓ stringr 1.4.0
✓ readr 1.4.0 ✓ forcats 0.5.1
── Conflicts ───────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
library(skimr)
Registered S3 methods overwritten by 'htmltools':
method from
print.html tools:rstudio
print.shiny.tag tools:rstudio
print.shiny.tag.list tools:rstudio
library(brms)
Loading required package: Rcpp
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Loading 'brms' package (version 2.15.7). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').
Attaching package: ‘brms’
The following object is masked from ‘package:stats’:
ar
anther <- read_csv("../input/Antherelongation2017.csv") %>%
rename_with(tolower) %>%
complete(day, floret, treatment, time)
── Column specification ───────────────────────────────────────────────────────────────────────────────────
cols(
Day = col_double(),
floret = col_double(),
Treatment = col_character(),
Time = col_time(format = ""),
`Length (mm)` = col_double()
)
anther <- anther %>% rename(trt=treatment, length=`length (mm)`) %>%
mutate(hour=as.numeric(time-min(time)) / 3600, # elapsed hours, first hour = 0
id=str_c(floret,day,trt,sep="-")) # unique id for each floret
anther %>% mutate(trt=as.factor(trt)) %>% skim()
── Data Summary ────────────────────────
Values
Name Piped data
Number of rows 1224
Number of columns 7
_______________________
Column type frequency:
character 1
difftime 1
factor 1
numeric 4
________________________
Group variables None
── Variable type: character ───────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate min max empty n_unique whitespace
1 id 0 1 8 12 0 72 0
── Variable type: difftime ────────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate min max median n_unique
1 time 0 1 19800 secs 34200 secs 27000 secs 17
── Variable type: factor ──────────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate ordered n_unique top_counts
1 trt 0 1 FALSE 3 Eas: 408, Wes: 408, Wes: 408
── Variable type: numeric ─────────────────────────────────────────────────────────────────────────────────
skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
1 day 0 1 5.12 2.71 1 2.75 5.5 7.25 9 ▇▃▃▇▇
2 floret 0 1 2 0.817 1 1 2 3 3 ▇▁▇▁▇
3 length 186 0.848 9.36 1.48 6.18 8.20 9.31 10.4 14.0 ▃▇▇▃▁
4 hour 0 1 2 1.23 0 1 2 3 4 ▇▆▆▆▇
unique(anther$trt)
[1] "East" "West" "WestHeat"
Plot to see what we have:
anther %>%
ggplot(aes(x=time, y=length, color=trt, shape=as.character(floret))) +
geom_point() +
facet_wrap(~day)
Warning: Removed 186 rows containing missing values (geom_point).
Fix height to max height. This is necessary because our growth model can’t account for shrinkage after reaching max height. Fill in missing end-of-series values, but note that they are censored.
anther <- anther %>%
mutate(censored = ifelse(is.na(length), "right", "none")) %>%
group_by(id) %>%
arrange(hour) %>%
mutate(length2=ifelse((length < max(length, na.rm = TRUE) | is.na(length)) & row_number() > which.max(length),
max(length, na.rm = TRUE), length)) %>%
ungroup() %>%
filter(!is.na(length2))
Plot smoothed data
anther %>%
ggplot(aes(x=time, y=length2, color=trt)) +
geom_smooth() +
facet_wrap(~day)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
And plot averaged over days:
anther %>%
ggplot(aes(x=time, y=length2, color=trt)) +
geom_smooth()
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
To help figure out reasonable priors
weibull <- function(Lmin, Lmax, k, delta, hour) {
y <- Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta)
tibble(hour=hour, y=y)
}
set up tibble with some priors
x <- seq(0,1,.01)
plot(x, dexp(x,1), type = "l")
k only
tibble(
Lmin=7.5,#rnorm(50, 7.5, 1),
Lmax=11.5,#rnorm(50, 11.5, 1),
k=rexp(50,1),
delta=1) %>% #rnorm(50,1,.5)) %>%
mutate(delta=ifelse(delta < 0, 0, delta)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
delta only
tibble(
Lmin=7.5,#rnorm(50, 7.5, 1),
Lmax=11.5,#rnorm(50, 11.5, 1),
k=.5,
delta=rexp(50,.5)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
delta and k
tibble(
Lmin=7.5,#rnorm(50, 7.5, 1),
Lmax=11.5,#rnorm(50, 11.5, 1),
k=rexp(100, 1),
delta=rexp(100,.5)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
all together
tibble(
Lmin=rnorm(100, 7.5, .25),
Lmax=rnorm(100, 11.5, .25),
k=rexp(100, 1),
delta=rexp(100, 0.5)) %>%
mutate(prediction=pmap(list(Lmin, Lmax, k, delta), weibull, seq(0,4.5,.1) ),
sample=row_number()) %>%
unnest(prediction) %>%
ggplot(aes(x=hour,y=y, group=sample)) +
geom_line(alpha=.2)
what priors are available to be set?
get_prior(bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin ~ 1,
delta + k ~ trt + (1|id),
nl=TRUE ),
data=anther)
prior class coef group resp dpar nlpar bound source
student_t(3, 0, 2.5) sigma default
(flat) b delta default
(flat) b Intercept delta (vectorized)
(flat) b trtWest delta (vectorized)
(flat) b trtWestHeat delta (vectorized)
student_t(3, 0, 2.5) sd delta default
student_t(3, 0, 2.5) sd id delta (vectorized)
student_t(3, 0, 2.5) sd Intercept id delta (vectorized)
(flat) b k default
(flat) b Intercept k (vectorized)
(flat) b trtWest k (vectorized)
(flat) b trtWestHeat k (vectorized)
student_t(3, 0, 2.5) sd k default
student_t(3, 0, 2.5) sd id k (vectorized)
student_t(3, 0, 2.5) sd Intercept id k (vectorized)
(flat) b Lmax default
(flat) b Intercept Lmax (vectorized)
(flat) b Lmin default
(flat) b Intercept Lmin (vectorized)
prior1 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k", lb=0),
prior(exponential(.5),nlpar="delta", lb=0) # 1 is exponential model
)
fit1c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin + k + delta ~ 1, #minimal model, parameters do not vary for treatment or random effects
nl=TRUE),
data=anther,
prior=prior1,
sample_prior = "yes",
cores=4,
chains=4)
Compiling Stan program...
Start sampling
starting worker pid=76107 on localhost:11119 at 12:58:08.890
starting worker pid=76121 on localhost:11119 at 12:58:09.129
starting worker pid=76135 on localhost:11119 at 12:58:09.348
starting worker pid=76149 on localhost:11119 at 12:58:09.601
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.001087 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.87 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2:
Chain 2: Gradient evaluation took 0.000662 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 6.62 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000723 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 7.23 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4:
Chain 4: Gradient evaluation took 0.00086 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 8.6 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 17.694 seconds (Warm-up)
Chain 1: 15.986 seconds (Sampling)
Chain 1: 33.68 seconds (Total)
Chain 1:
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 17.407 seconds (Warm-up)
Chain 2: 16.899 seconds (Sampling)
Chain 2: 34.306 seconds (Total)
Chain 2:
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 17.649 seconds (Warm-up)
Chain 3: 17.756 seconds (Sampling)
Chain 3: 35.405 seconds (Total)
Chain 3:
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 23.21 seconds (Warm-up)
Chain 4: 13.42 seconds (Sampling)
Chain 4: 36.63 seconds (Total)
Chain 4:
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4216131 225.2 6937521 370.6 NA 4969379 265.4
Vcells 8915431 68.1 14964384 114.2 16384 12403654 94.7
fit1c <- add_criterion(fit1c, "loo")
summary(fit1c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
k ~ 1
delta ~ 1
Data: anther (Number of observations: 1224)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 12.51 0.15 12.24 12.82 1.00 2167 2304
Lmin_Intercept 7.94 0.07 7.79 8.08 1.00 2539 2404
k_Intercept 0.39 0.01 0.37 0.42 1.00 2181 2271
delta_Intercept 2.29 0.14 2.03 2.59 1.00 2254 2127
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 1.12 0.02 1.07 1.17 1.00 3041 2272
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4258160 227.5 6937521 370.6 NA 6829137 364.8
Vcells 9011615 68.8 50512828 385.4 16384 63141016 481.8
plot(fit1c)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4388888 234.4 6937521 370.6 NA 6937521 370.6
Vcells 9851967 75.2 40410263 308.4 16384 63141016 481.8
pairs(fit1c)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4572908 244.3 6937521 370.6 NA 6937521 370.6
Vcells 10783056 82.3 40410263 308.4 16384 63141016 481.8
pred1c <- predict(fit1c)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4286391 229.0 6937521 370.6 NA 6937521 370.6
Vcells 9064549 69.2 60381579 460.7 16384 75476557 575.9
pred1c <- cbind(pred1c, anther)
pred1c %>% ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate),color="yellow")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pp_check(fit1c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
prior2 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k", lb=0),
prior(exponential(0.5),nlpar="delta"), # 1 is exponential model
prior(normal(0,0.5), nlpar="delta", coef="trtWest"),
prior(normal(0,0.5), nlpar="delta", coef="trtWestHeat")
)
fit2c <- brm( # ignore the warning about lower bounds. If we set it for the intercept, it impacts the other delta coefficients and we don't want that.
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin + k ~ 1,
delta ~ trt + (1|id), #
nl=TRUE ),
data=anther,
prior=prior2,
sample_prior = "yes",
cores=4,
chains=4,
iter=5000 #,
#control=list(adapt_delta=0.9)
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=76306 on localhost:11119 at 13:00:33.794
starting worker pid=76320 on localhost:11119 at 13:00:34.073
starting worker pid=76334 on localhost:11119 at 13:00:34.333
starting worker pid=76348 on localhost:11119 at 13:00:34.589
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1:
Chain 1: Gradient evaluation took 0.000766 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.66 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -1.84619, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -1.51963, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -0.609171, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -1.11426, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -1.24104, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -0.099825, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -1.49976, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2: Rejecting initial value:
Chain 2: Error evaluating the log probability at the initial value.
Chain 2: Exception: exponential_lpdf: Random variable is -1.90058, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 2:
Chain 2: Gradient evaluation took 0.00073 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 7.3 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.002524 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 25.24 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: exponential_lpdf: Random variable is -0.26906, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: exponential_lpdf: Random variable is -1.68216, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 4: Rejecting initial value:
Chain 4: Error evaluating the log probability at the initial value.
Chain 4: Exception: exponential_lpdf: Random variable is -1.23181, but must be nonnegative! (in 'anon_model', line 78, column 2 to column 47)
Chain 4:
Chain 4: Gradient evaluation took 0.001643 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 16.43 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 145.761 seconds (Warm-up)
Chain 4: 119.81 seconds (Sampling)
Chain 4: 265.571 seconds (Total)
Chain 4:
Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 143.16 seconds (Warm-up)
Chain 1: 141.503 seconds (Sampling)
Chain 1: 284.663 seconds (Total)
Chain 1:
Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 144.213 seconds (Warm-up)
Chain 2: 152.251 seconds (Sampling)
Chain 2: 296.464 seconds (Total)
Chain 2:
Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 214.18 seconds (Warm-up)
Chain 3: 97.019 seconds (Sampling)
Chain 3: 311.199 seconds (Total)
Chain 3:
fit2c <- add_criterion(fit2c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4498318 240.3 6937521 370.6 NA 6937521 370.6
Vcells 11825816 90.3 126399944 964.4 16384 157999363 1205.5
summary(fit2c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
k ~ 1
delta ~ trt + (1 | id)
Data: anther (Number of observations: 1224)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
Group-Level Effects:
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 0.97 0.11 0.78 1.20 1.00 1845 3119
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 15.51 0.14 15.22 15.80 1.00 5315 5822
Lmin_Intercept 7.80 0.05 7.70 7.90 1.00 7920 7043
k_Intercept 0.24 0.00 0.23 0.25 1.00 5199 6148
delta_Intercept 1.75 0.19 1.39 2.14 1.00 1240 2378
delta_trtWest 0.19 0.25 -0.29 0.68 1.00 1277 2338
delta_trtWestHeat 0.16 0.24 -0.33 0.63 1.00 1353 2573
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.77 0.02 0.74 0.81 1.00 7473 7308
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4498409 240.3 6937521 370.6 NA 6937521 370.6
Vcells 11826074 90.3 101119956 771.5 16384 157999363 1205.5
plot(fit2c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4689628 250.5 6937521 370.6 NA 6937521 370.6
Vcells 14290147 109.1 80895965 617.2 16384 157999363 1205.5
plot(conditional_effects(fit2c), ask=FALSE)
pp_check(fit2c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred2c <- predict(fit2c)
pred2c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred2c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales = "free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred2c %>%
cbind(anther) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c)
elpd_diff se_diff
fit2c 0.0 0.0
fit1c -376.6 22.1
Fit2c is preferred, but predictions are not that good relative to observed
hypothesis(fit2c, c("delta_trtWest = 0",
"delta_trtWestHeat = 0",
"delta_trtWest = delta_trtWestHeat"))
Hypothesis Tests for class b:
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.
inits3 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=rexp(1, .5),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(72,1362))
)
}
prior3 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.02), nlpar="k", coef="trtWest"),
prior(normal(0,0.02), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5),nlpar="delta")
)
fit3c <- brm( # ignore the warning about k.
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin + delta ~ 1,
k ~ trt + (1|id), #
nl=TRUE ),
data=anther,
prior=prior3,
sample_prior = "yes",
inits = inits3,
cores=4,
chains=4,
iter=5000,
control=list(adapt_delta=0.99)
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=76542 on localhost:11119 at 13:07:47.250
starting worker pid=76556 on localhost:11119 at 13:07:47.528
starting worker pid=76570 on localhost:11119 at 13:07:47.766
starting worker pid=76584 on localhost:11119 at 13:07:48.013
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000538 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.38 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.0036 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 36 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.001676 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 16.76 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.000881 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 8.81 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 262.764 seconds (Warm-up)
Chain 3: 191.639 seconds (Sampling)
Chain 3: 454.403 seconds (Total)
Chain 3:
Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 276.656 seconds (Warm-up)
Chain 1: 187.693 seconds (Sampling)
Chain 1: 464.349 seconds (Total)
Chain 1:
Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 280.082 seconds (Warm-up)
Chain 2: 190.691 seconds (Sampling)
Chain 2: 470.773 seconds (Total)
Chain 2:
Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 273.483 seconds (Warm-up)
Chain 4: 220.496 seconds (Sampling)
Chain 4: 493.979 seconds (Total)
Chain 4:
fit3c <- add_criterion(fit3c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4720647 252.2 6937521 370.6 NA 6937521 370.6
Vcells 15303922 116.8 138496461 1056.7 16384 173119397 1320.8
summary(fit3c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
delta ~ 1
k ~ trt + (1 | id)
Data: anther (Number of observations: 1224)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
Group-Level Effects:
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(k_Intercept) 0.13 0.01 0.11 0.16 1.00 1060 2298
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 13.15 0.12 12.93 13.40 1.00 5403 5601
Lmin_Intercept 7.78 0.05 7.68 7.88 1.00 6050 6812
delta_Intercept 1.95 0.08 1.80 2.12 1.00 4847 6181
k_Intercept 0.37 0.02 0.34 0.41 1.01 695 1735
k_trtWest -0.01 0.02 -0.04 0.03 1.00 1557 3077
k_trtWestHeat -0.00 0.02 -0.04 0.03 1.00 1786 3598
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.72 0.02 0.69 0.75 1.00 8514 7420
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4720712 252.2 6937521 370.6 NA 6937521 370.6
Vcells 15304127 116.8 110797169 845.4 16384 173119397 1320.8
pairs(fit3c, pars="^b_")
plot(fit3c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4913619 262.5 9099172 486.0 NA 6937521 370.6
Vcells 17821208 136.0 70910189 541.1 16384 173119397 1320.8
plot(conditional_effects(fit3c), ask=FALSE)
pp_check(fit3c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred3c <- predict(fit3c)
pred3c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred3c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales = "free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred3c %>%
cbind(anther) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c, fit3c)
elpd_diff se_diff
fit3c 0.0 0.0
fit2c -81.6 19.2
fit1c -458.2 24.0
inits4 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=c(rexp(1, .5), 0, 0),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(72,1)),
z_2=array(0, dim=c(72,1))
)
}
prior4 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.02), nlpar="k", coef="trtWest"),
prior(normal(0,0.02), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5), nlpar="delta"), # 1 is exponential model
prior(normal(0,0.25), nlpar="delta", coef="trtWest"),
prior(normal(0,0.25), nlpar="delta", coef="trtWestHeat")
)
fit_4c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin ~ 1,
delta + k ~ trt + (1|id), #
nl=TRUE ),
data=anther,
prior=prior4,
inits=inits4,
cores=4,
chains=4,
iter=5000, control = list(adapt_delta=0.95)
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=76837 on localhost:11119 at 13:18:05.011
starting worker pid=76851 on localhost:11119 at 13:18:05.276
starting worker pid=76865 on localhost:11119 at 13:18:05.520
starting worker pid=76879 on localhost:11119 at 13:18:05.768
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000894 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 8.94 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2:
Chain 2: Gradient evaluation took 0.000906 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 9.06 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.000832 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 8.32 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4: Log probability evaluates to log(0), i.e. negative infinity.
Chain 4: Stan can't start sampling from this initial value.
Chain 4:
Chain 4: Gradient evaluation took 0.012245 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 122.45 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 330.179 seconds (Warm-up)
Chain 1: 235.254 seconds (Sampling)
Chain 1: 565.433 seconds (Total)
Chain 1:
Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 332.767 seconds (Warm-up)
Chain 2: 234.393 seconds (Sampling)
Chain 2: 567.16 seconds (Total)
Chain 2:
Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 331.984 seconds (Warm-up)
Chain 3: 236.42 seconds (Sampling)
Chain 3: 568.404 seconds (Total)
Chain 3:
Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 333.266 seconds (Warm-up)
Chain 4: 235.153 seconds (Sampling)
Chain 4: 568.419 seconds (Total)
Chain 4:
fit_4c <- add_criterion(fit_4c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4924736 263.1 9099172 486.0 NA 9099172 486.0
Vcells 19376727 147.9 121404308 926.3 16384 189693904 1447.3
summary(fit_4c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ 1
Lmin ~ 1
delta ~ trt + (1 | id)
k ~ trt + (1 | id)
Data: anther (Number of observations: 1224)
Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
total post-warmup samples = 10000
Group-Level Effects:
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 0.58 0.09 0.42 0.76 1.00 3963 6578
sd(k_Intercept) 0.12 0.01 0.10 0.14 1.00 2479 4528
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 13.40 0.13 13.15 13.67 1.00 9347 8262
Lmin_Intercept 7.76 0.05 7.68 7.85 1.00 16916 7639
delta_Intercept 1.79 0.14 1.53 2.08 1.00 5701 6731
delta_trtWest 0.29 0.16 -0.03 0.60 1.00 6782 7010
delta_trtWestHeat 0.35 0.16 0.03 0.64 1.00 6536 7190
k_Intercept 0.36 0.02 0.32 0.39 1.00 1951 4043
k_trtWest -0.00 0.02 -0.04 0.03 1.00 3971 6113
k_trtWestHeat -0.00 0.02 -0.03 0.03 1.00 3906 6122
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.67 0.02 0.64 0.70 1.00 13486 6615
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 4924796 263.1 9099172 486 NA 9099172 486.0
Vcells 19376968 147.9 97123447 741 16384 189693904 1447.3
pairs(fit_4c, pars = "^b_")
plot(fit_4c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5099113 272.4 9099172 486.0 NA 9099172 486.0
Vcells 22025898 168.1 77698758 592.8 16384 189693904 1447.3
plot(conditional_effects(fit_4c), ask=FALSE)
pp_check(fit_4c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred_4c <- predict(fit_4c)
pred_4c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_4c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales="free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_4c %>%
cbind(anther) %>%
group_by(trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_line() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2)
`summarise()` has grouped output by 'trt'. You can override using the `.groups` argument.
pred_4c %>%
cbind(anther) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c, fit3c, fit_4c)
elpd_diff se_diff
fit_4c 0.0 0.0
fit3c -68.5 9.8
fit2c -150.1 19.2
fit1c -526.7 26.1
inits5 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=c(rexp(1, .5), 0, 0),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(8,1)),
z_2=array(0, dim=c(72,1)),
z_3=array(0, dim=c(72,1))
)
}
prior5 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.02), nlpar="k", coef="trtWest"),
prior(normal(0,0.02), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5), nlpar="delta"), # 1 is exponential model
prior(normal(0,0.25), nlpar="delta", coef="trtWest"),
prior(normal(0,0.25), nlpar="delta", coef="trtWestHeat")
)
fit_5c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmin ~ 1,
Lmax ~ (1|day),
delta + k ~ trt + (1|id), #
nl=TRUE ),
data=anther,
prior=prior5,
sample_prior = "yes",
inits=inits5,
cores=4,
chains=4,
iter=2000
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=77180 on localhost:11119 at 13:30:40.708
starting worker pid=77195 on localhost:11119 at 13:30:41.020
starting worker pid=77209 on localhost:11119 at 13:30:41.363
starting worker pid=77224 on localhost:11119 at 13:30:41.846
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.001087 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.87 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 0.001937 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 19.37 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3: Log probability evaluates to log(0), i.e. negative infinity.
Chain 3: Stan can't start sampling from this initial value.
Chain 3:
Chain 3: Gradient evaluation took 0.065723 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 657.23 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.00087 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 8.7 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 254.161 seconds (Warm-up)
Chain 4: 97.843 seconds (Sampling)
Chain 4: 352.004 seconds (Total)
Chain 4:
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 266.991 seconds (Warm-up)
Chain 1: 98.1 seconds (Sampling)
Chain 1: 365.091 seconds (Total)
Chain 1:
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 259.674 seconds (Warm-up)
Chain 3: 104.776 seconds (Sampling)
Chain 3: 364.45 seconds (Total)
Chain 3:
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 253.592 seconds (Warm-up)
Chain 2: 126.031 seconds (Sampling)
Chain 2: 379.623 seconds (Total)
Chain 2:
fit_5c <- add_criterion(fit_5c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5130163 274.0 9099172 486.0 NA 9099172 486.0
Vcells 21717648 165.7 119569895 912.3 16384 189693904 1447.3
summary(fit_5c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmin ~ 1
Lmax ~ (1 | day)
delta ~ trt + (1 | id)
k ~ trt + (1 | id)
Data: anther (Number of observations: 1224)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~day (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Lmax_Intercept) 2.75 0.81 1.55 4.65 1.00 729 1196
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 0.67 0.11 0.48 0.91 1.00 875 1567
sd(k_Intercept) 0.04 0.01 0.03 0.05 1.00 987 2173
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmin_Intercept 7.73 0.05 7.63 7.82 1.00 3804 3102
Lmax_Intercept 11.65 0.26 11.14 12.14 1.00 1669 1936
delta_Intercept 1.76 0.16 1.47 2.08 1.00 1018 1553
delta_trtWest 0.23 0.17 -0.09 0.56 1.00 1182 2008
delta_trtWestHeat 0.30 0.17 -0.04 0.62 1.00 1244 1628
k_Intercept 0.33 0.02 0.28 0.36 1.00 893 1650
k_trtWest -0.00 0.01 -0.03 0.02 1.00 1747 2496
k_trtWestHeat 0.00 0.01 -0.02 0.03 1.00 1594 2321
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.66 0.02 0.64 0.69 1.00 3332 2873
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5130222 274.0 9099172 486.0 NA 9099172 486.0
Vcells 21717889 165.7 95655916 729.8 16384 189693904 1447.3
plot(fit_5c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5316379 284.0 9099172 486.0 NA 9099172 486.0
Vcells 23223721 177.2 76524733 583.9 16384 189693904 1447.3
plot(conditional_effects(fit_5c), ask=FALSE)
pp_check(fit_5c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred_5c <- predict(fit_5c)
pred_5c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_5c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_5c %>%
cbind(anther) %>%
group_by(trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_line() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2)
`summarise()` has grouped output by 'trt'. You can override using the `.groups` argument.
pred_5c %>%
cbind(anther) %>%
group_by(day, trt, hour) %>%
summarize(length2 = mean(length2),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
loo_compare(fit1c, fit2c, fit3c, fit_4c, fit_5c)
elpd_diff se_diff
fit_5c 0.0 0.0
fit_4c -21.9 6.7
fit3c -90.4 13.0
fit2c -172.0 21.1
fit1c -548.6 28.0
inits6 <- function() { list(
b_Lmax=rnorm(1, 11.5, .25),
b_Lmin=rnorm(1, 7.5, .25),
b_delta=c(rexp(1, .5), 0, 0),
b_k=c(rexp(1, 1), 0, 0),
z_1=array(0, dim=c(8,1)),
z_2=array(0, dim=c(8,1)),
z_3=array(0, dim=c(72,1)),
z_4=array(0, dim=c(72,1))
)
}
prior6 <- c(prior(normal(7.5,.25), nlpar="Lmin"),
prior(normal(11.5, .25),nlpar="Lmax"),
prior(exponential(1),nlpar="k"),
prior(normal(0,0.04), nlpar="k", coef="trtWest"),
prior(normal(0,0.04), nlpar="k", coef="trtWestHeat"),
prior(exponential(0.5), nlpar="delta"), # 1 is exponential model
prior(normal(0,0.25), nlpar="delta", coef="trtWest"),
prior(normal(0,0.25), nlpar="delta", coef="trtWestHeat")
)
fit_6c <- brm(
formula=bf(
length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k*hour)^delta), #Weibull model
Lmax + Lmin ~ (1|day),
delta + k ~ trt + (1|id), #
nl=TRUE ),
data=anther,
prior=prior6,
sample_prior = "yes",
inits=inits6,
cores=4,
chains=4,
iter=2000
)
Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
Warning occurred for prior
b_k ~ exponential(1)
b_delta ~ exponential(0.5)
Compiling Stan program...
Start sampling
starting worker pid=77463 on localhost:11119 at 13:39:02.227
starting worker pid=77477 on localhost:11119 at 13:39:02.504
starting worker pid=77491 on localhost:11119 at 13:39:02.754
starting worker pid=77505 on localhost:11119 at 13:39:03.012
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can't start sampling from this initial value.
Chain 1:
Chain 1: Gradient evaluation took 0.001143 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.43 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2: Log probability evaluates to log(0), i.e. negative infinity.
Chain 2: Stan can't start sampling from this initial value.
Chain 2:
Chain 2: Gradient evaluation took 0.001321 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 13.21 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 0.001497 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 14.97 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 0.001221 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 12.21 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 224.159 seconds (Warm-up)
Chain 4: 84.721 seconds (Sampling)
Chain 4: 308.88 seconds (Total)
Chain 4:
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 230.212 seconds (Warm-up)
Chain 3: 83.288 seconds (Sampling)
Chain 3: 313.5 seconds (Total)
Chain 3:
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 241.278 seconds (Warm-up)
Chain 1: 79.375 seconds (Sampling)
Chain 1: 320.653 seconds (Total)
Chain 1:
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 242.74 seconds (Warm-up)
Chain 2: 78.142 seconds (Sampling)
Chain 2: 320.882 seconds (Total)
Chain 2:
fit_6c <- add_criterion(fit_6c, "loo")
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5334636 285.0 9099172 486.0 NA 9099172 486.0
Vcells 24170397 184.5 76615362 584.6 16384 189693904 1447.3
summary(fit_6c)
Family: gaussian
Links: mu = identity; sigma = identity
Formula: length2 | cens(censored) ~ Lmax - (Lmax - Lmin) * exp(-(k * hour)^delta)
Lmax ~ (1 | day)
Lmin ~ (1 | day)
delta ~ trt + (1 | id)
k ~ trt + (1 | id)
Data: anther (Number of observations: 1224)
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup samples = 4000
Group-Level Effects:
~day (Number of levels: 8)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Lmax_Intercept) 2.47 0.72 1.39 4.16 1.00 1024 1532
sd(Lmin_Intercept) 0.61 0.22 0.33 1.14 1.00 970 1719
~id (Number of levels: 72)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(delta_Intercept) 0.49 0.10 0.32 0.71 1.00 949 1722
sd(k_Intercept) 0.04 0.01 0.03 0.05 1.00 1013 1946
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Lmax_Intercept 11.66 0.24 11.17 12.15 1.00 1903 2592
Lmin_Intercept 7.71 0.17 7.36 8.02 1.00 1069 2003
delta_Intercept 1.76 0.15 1.48 2.06 1.00 1206 1981
delta_trtWest 0.35 0.15 0.06 0.64 1.00 1950 2516
delta_trtWestHeat 0.41 0.15 0.10 0.69 1.00 1799 2631
k_Intercept 0.34 0.02 0.29 0.38 1.00 948 1398
k_trtWest -0.00 0.01 -0.03 0.02 1.00 1415 2297
k_trtWestHeat 0.00 0.01 -0.03 0.03 1.00 1443 1947
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.64 0.02 0.61 0.67 1.00 2616 2625
Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5334708 285.0 9099172 486.0 NA 9099172 486.0
Vcells 24170671 184.5 76615362 584.6 16384 189693904 1447.3
plot(fit_6c, ask = FALSE)
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5538949 295.9 9099172 486.0 NA 9099172 486.0
Vcells 25827427 197.1 76615362 584.6 16384 189693904 1447.3
plot(conditional_effects(fit_6c), ask=FALSE)
pp_check(fit_6c)
Using 10 posterior samples for ppc type 'dens_overlay' by default.
Warning: Censored responses are not shown in 'pp_check'.
pred_6c <- predict(fit_6c)
pred_6c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length2, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1)
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_6c %>%
cbind(anther) %>%
ggplot(aes(x=hour, y=length, color=trt)) +
geom_smooth(span=1) +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_smooth(aes(y=Estimate), lty=2, span=1) +
facet_wrap(~day, scales = "free_y")
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
Warning: Removed 186 rows containing non-finite values (stat_smooth).
`geom_smooth()` using method = 'loess' and formula 'y ~ x'
pred_6c %>%
cbind(anther) %>%
group_by(day, trt, hour) %>%
summarize(length = mean(length),
Estimate = mean(Estimate)) %>%
ggplot(aes(x=hour, y=length, color=trt)) +
geom_point() +
#geom_ribbon(aes(ymin=Q2.5, ymax=Q97.5), color="grey90", fill="black", alpha=.1) +
geom_line(aes(y=Estimate), lty=2) +
facet_wrap(~day, scales = "free_y", ncol = 4) +
theme(legend.position = "top")
`summarise()` has grouped output by 'day', 'trt'. You can override using the `.groups` argument.
Warning: Removed 62 rows containing missing values (geom_point).
ggsave("../output/anther_fit_plot.pdf", width=6.5, height = 4)
Warning: Removed 62 rows containing missing values (geom_point).
loo_compare(fit1c, fit2c, fit3c, fit_4c, fit_5c, fit_6c)
elpd_diff se_diff
fit_6c 0.0 0.0
fit_5c -41.2 13.3
fit_4c -63.1 14.7
fit3c -131.6 17.1
fit2c -213.2 21.4
fit1c -589.8 27.2
Get the posterior
post_6c <- posterior_samples(fit_6c, pars="b.*")
posterior for inflection point
# inflection point equation
infl.pt <- function(delta, k) {
(1/k) *
((delta -1) / delta)^(1/delta)
}
post_6c$infl_intercept <- infl.pt(post_6c$b_delta_Intercept, post_6c$b_k_Intercept)
post_6c$infl_trtWest <-
infl.pt(
post_6c$b_delta_Intercept + post_6c$b_delta_trtWest,
post_6c$b_k_Intercept + post_6c$b_k_trtWest)
post_6c$infl_trtWestHeat <-
infl.pt(
post_6c$b_delta_Intercept+post_6c$b_delta_trtWestHeat,
post_6c$b_k_Intercept + post_6c$b_k_trtWestHeat)
apply(post_6c, 2, function(x) mean(x))
b_Lmax_Intercept b_Lmin_Intercept b_delta_Intercept b_delta_trtWest
11.6603550955 7.7107106754 1.7572125610 0.3480487032
b_delta_trtWestHeat b_k_Intercept b_k_trtWest b_k_trtWestHeat
0.4085974213 0.3366791988 -0.0043306062 0.0021603672
prior_b_Lmax prior_b_Lmin prior_b_delta_Intercept prior_b_delta_trtWest
11.4986398949 7.4999693078 2.0273121585 -0.0011335972
prior_b_delta_trtWestHeat prior_b_k_Intercept prior_b_k_trtWest prior_b_k_trtWestHeat
-0.0038707272 1.0067453076 0.0003744252 0.0011928831
infl_intercept infl_trtWest infl_trtWestHeat
1.8181198665 2.2064981348 2.2098722873
posterior probability that inflection point for trtWestHeat is greater than inflection point for East
sum( (post_6c$infl_trtWestHeat > post_6c$infl_intercept )) / nrow(post_6c)
[1] 0.99
posterior probability that inflection point for trtWest is greater than inflection point for East
sum( (post_6c$infl_trtWest > post_6c$infl_intercept )) / nrow(post_6c)
[1] 0.991
posterior probability that inflection point for trtWest is greater than inflection point for trtWestHeat
sum((post_6c$infl_trtWest > post_6c$infl_trtWestHeat)) / nrow(post_6c)
[1] 0.49175
posterior probability that k for trtWestHeat is less than East
sum( (post_6c$b_k_trtWestHeat < 0 )) / nrow(post_6c)
[1] 0.43225
sum( (post_6c$b_k_trtWestHeat + post_6c$b_k_Intercept < post_6c$b_k_Intercept )) / nrow(post_6c)
[1] 0.43225
posterior probability that k for trtWest is less than East
sum( (post_6c$b_k_trtWest < 0 )) / nrow(post_6c)
[1] 0.62575
sum( (post_6c$b_k_trtWest + post_6c$b_k_Intercept < post_6c$b_k_Intercept )) / nrow(post_6c)
[1] 0.62575
posterior probability that k for trtWest is less than trtWestHeat
sum( (post_6c$b_k_trtWest < post_6c$b_k_trtWestHeat )) / nrow(post_6c)
[1] 0.701
k
post_6c %>%
as.data.frame() %>%
select(starts_with("b_k")) %>%
rename_with(~str_remove(.,"b_")) %>%
rename_with(~str_remove(.,"trt")) %>%
rename(k_East=k_Intercept) %>%
mutate(k_West=k_West+k_East,
k_WestHeat=k_WestHeat+k_East) %>%
summarize(across(.fns = list(mean = mean,
lower = ~ rethinking::HPDI(., 0.95)[1],
upper = ~ rethinking::HPDI(., 0.95)[2])
)) %>%
pivot_longer(cols=everything(), names_to=c("trt", "parameter"), names_pattern = "k_(.*)_(.*)") %>%
pivot_wider(names_from = "parameter", values_from = "value") %>%
ggplot(aes(x=trt, y=mean, ymin=lower, ymax=upper)) +
geom_col(fill="skyblue") +
geom_errorbar(width=.5) +
ylab("k") +
xlab("Condition") +
theme_bw()
ggsave("../output/anther_k_plot.pdf", height=3, width = 3)
inflection
post_6c %>%
as.data.frame() %>%
select(starts_with("infl")) %>%
rename_with(~str_remove(.,"trt")) %>%
rename_with(~str_remove(.,"infl_")) %>%
rename(East=intercept) %>%
summarize(across(.fns = list(mean = mean,
lower = ~ rethinking::HPDI(., 0.95)[1],
upper = ~ rethinking::HPDI(., 0.95)[2])
)) %>%
pivot_longer(cols=everything(), names_to=c("trt", "parameter"), names_pattern = "(.*)_(.*)") %>%
pivot_wider(names_from = "parameter", values_from = "value") %>%
ggplot(aes(x=trt, y=mean, ymin=lower, ymax=upper)) +
geom_col(fill="red") +
geom_errorbar(width=.5) +
ylab("Inflection Point") +
xlab("Condition") +
theme_bw()
ggsave("../output/anther_inflection_plot.pdf", height=3, width = 3)